Advances in flow and mass cytometry now enable the simultaneous measurement of more than 50 parameters per cell, transforming it into a high-dimensional technology. Together with the high throughput of analyzing millions of cells per sample, the data requires sophisticated analysis techniques. Conventional gating, until now the gold standard for cytometry data analysis, faces challenges in scaling for the application of high-dimensional cytometry data. Issues with conventional gating strategies include inefficiency, subjectivity, and a significant risk of missing unknown or minor cell populations. Hence, unbiased and computational methods for high-dimensional cytometry data analysis are required.
Several computational tools for cytometry data analysis are available, primarily designed for the R programming language. Examples include CATALYST / CyTOF workflow, FlowSOM, flowCore, flowAI, and diffcyt. These tools address specific aspects of high-dimensional cytometry data analysis, including handling flow cytometry standard (FCS) files, preprocessing, dimensionality reduction, clustering based on self-organizing maps, and differential abundance analysis. However, interoperability between these tools is often limited and inefficient, and methods for additional in-depth data analysis are lacking.
Whereas high-dimensional cytometry data analysis methods are less established, tools for analyzing single-cell RNA sequencing (scRNAseq) data are highly developed and offer in-depth analysis approaches. The most commonly used framework for scRNAseq data analysis is Seurat, which covers every aspect of data analysis from preprocessing to a range of downstream analyses, including dimensionality reduction and clustering. A notable advantage of scRNAseq tools is their interoperability. Most third-party tools that offer additional analysis methods, such as batch correction or trajectory analyses, integrate seamlessly into the Seurat framework. However, these scRNAseq analysis tools are not readily accessible for the analysis of cytometry data.
Our novel and user-friendly high-dimensional cytometry analysis R framework, Seumetry, utilizes a broad range of cytometry-specific tools for data handling and integrates into Seurat, providing access to additional scRNAseq analysis tools. The workflow comprises loading of FCS files, preprocessing (compensation and transformation), quality control, and downstream analysis in Seurat, such as dimensionality reduction, clustering, and more. Our advanced quality control algorithms cover the detection of outliers and antibody aggregate artifacts. Antibody aggregates may occur with highly complex cytometry panels, leading to false positive events. Adjusting laboratory protocols can mitigate these aggregates, but these adjustments are often impractical and a lack of rare samples may prohibit re-analysis. Our novel algorithm detects these aggregate artifacts and cleans the data. For downstream analysis, we implemented additional methods covering differential expression and abundance analysis, as well as visualization techniques. In summary, Seumetry combines all aspects of cytometry data analysis in a user-friendly R package, comprising loading, preprocessing, quality control, differential expression and abundance, visualizations, and downstream analyses through the Seurat framework.
To demonstrate the functionality and workflow of Seumetry, we employed a high-dimensional spectral flow cytometry dataset comprising more than 40 parameters of human intestinal immune cells. The dataset consists of 10 intestinal mucosal samples (epithelium = IEL, and lamina propria = LPL) from 5 adult human donors. The raw data was unmixed and pre-gated on single, live CD45-positive immune cells.
Load libraries
library(Seumetry)
library(ggplot2)
library(dplyr)
Set global ggplot2 options and define theme to make beautiful plots
# adjust fill and color globally using options
discrete_colors <- c("#5988B2", "#C9635E", "#67976B", "#C88F67", "#A583B0", "#CDB86A",
"#ABCCE2", "#7F007F", "#C8A5A3", "#B9E1B8", "#C5807B", "#E0CDB1", "#917C6F",
"#5988B2", "#D5E8AE", "#CC7E78", "#B4C9E3", "#AB8BAF", "#97C496", "#5F6E9D",
"#B3CC96", "#AF8B99", "#C88F67", "#C16D6B", "#8E8E8E", "#5E93C2", "#FF8FBF",
"#AFD192", "#8E8E8E", "#8BCF92", "#5F6E9D", "#8E6A8E", "#A88169", "#9D7BAE",
"#6D8D7E", "#6E6E6E", "#9D4949", "#F2E89C", "#A9FF93", "#93A9C5", "#B76CA8",
"#D4857B", "#4C0000", "#BDBA88", "#6E6E6E", "#FFDF7F", "#B2B2B2", "#366E5E",
"#C6C6C6", "#DDBE8F", "#C6C6C6", "#89603A", "#9C7A8D", "#7CCCBE", "#AFFF9F")
options(ggplot2.discrete.colour = discrete_colors)
options(ggplot2.discrete.fill = discrete_colors)
# define new gradient colors
gradient_colors <- colorRampPalette(c("steelblue", "slategray2", "white", "tan1",
"firebrick3"))(100)
# define a theme to match plots from other packages with Seumetry plots
set_theme <- list(theme_linedraw(), theme(aspect.ratio = 1, panel.grid.major = element_blank(),
panel.grid.minor = element_blank()))
Metadata and panel are data.frames required by Seumetry for loading FCS files, creating a Seurat object, and preprocess the data.
metadata <- read.csv("data/metadata.csv")
panel <- read.csv("data/panel.csv")
Metadata: a data.frame with at least 2 columns: “file_name” and “sample_id”. All additional columns are used as metadata columns and are added to the Seurat object.
head(metadata)
## file_name sample_id patient_id tissue layer condition
## 1 303_Ileum_IEL_CD45pos.fcs 303_Ileum_IEL 303 ileum IEL healthy
## 2 303_Ileum_LPL_CD45pos.fcs 303_Ileum_LPL 303 ileum LPL healthy
## 3 304_Ileum_IEL_CD45pos.fcs 304_Ileum_IEL 304 ileum IEL healthy
## 4 304_Ileum_LPL_CD45pos.fcs 304_Ileum_LPL 304 ileum LPL healthy
## 5 306_Ileum_IEL_CD45pos.fcs 306_Ileum_IEL 306 ileum IEL healthy
## 6 306_Ileum_LPL_CD45pos.fcs 306_Ileum_LPL 306 ileum LPL healthy
Panel: a data.frame with at least 2 columns: “fcs_colname” and “antigen”. The fcs_colname are the names of the channels in the FCS file. Antigen is the desired name for downstream analysis. Additionally, columns for transformation are required. See below (transformation) for details.
head(panel)
## fcs_colname antigen arcsinh_cofactor biexp_neg biexp_pos biexp_width
## 1 APC.A IL15RA 6000 0 4.0 -1000
## 2 APC.Cy7.A CD27 6000 0 4.5 -1000
## 3 APC.Fire.810.A CCR7 6000 0 4.0 -1000
## 4 APC.R700.A CD127 6000 0 4.5 -1000
## 5 Alexa.Fluor.647.A CD1C 6000 1 4.5 -1000
## 6 BB515.A CD141 6000 1 6.0 -1000
FCS files are loaded based on *.fcs files in the folder and on “file_name” column in metadata data.frame.
Each cell will receive a unique cell name based on sample_id provided in metadata.
fcs_fs <- create_flowset("data/fcs", metadata)
fcs_fs
## A flowSet with 10 experiments.
##
## column names(50): FSC.A FSC.H ... eFluor.506.A Time
It is highly recommended to save the “fcs_fs” object, as it can be used later to export FCS files based on specific cells.
saveRDS(fcs_fs, "fcs_fs.rds")
All FCS files will be merged and metadata will be added based on sample_id and metadata data.frame. Raw data is saved in Seurat assay “fcs”. Only Channels will be kept that are present in panel data.frame and channels will be renamed from “fcs_colname” to “antigen”.
All panel data are stored in seu@misc slot for easy access.
seu <- create_seurat(fcs_fs, panel, metadata)
seu
## An object of class Seurat
## 50 features across 1124603 samples within 2 assays
## Active assay: fcs (39 features, 0 variable features)
## 1 other assay present: unused
# plot cellnumbers per sample
plot_cellnumber(seu)
Compensation is based on spillover matrix that can be provided within the FCS file or as a data.frame.
Option 1) Use external spillover matrix directly: same compensation for all files.
Option 2) Use spillover matrix saved in FCS files: same compensation for all files.
Option 3) Use spillover matrix saved in FCS files: compensation matrix used from individual FCS files, thus can be different for each file.
The compensate_data function will use flowCore::compensate to compensate the raw values and write a new assay into the Seurat object called “comp”. The compensated data are saved in “counts” slot.
Warning: if multiple spillover matrices are present in FCS files, use the correct one!
Here, we use option 3 to compensate the data.
# Check different matrices using:
names(flowCore::spillover(fcs_fs[[1]]))
## [1] "SPILL" "spillover" "$SPILLOVER"
# In this case, spillover matrix in column 3 is the correct one.
seu <- compensate_data(fcs_fs, seu, fcs_matrix = 3)
# Check that compensation worked.
plot1 <- plot_cyto(seu, x = "IL15-2RB", y = "CD69", assay = "fcs", slot = "counts",
scale = "log") + ggtitle("Uncompensated")
plot2 <- plot_cyto(seu, x = "IL15-2RB", y = "CD69", assay = "comp", slot = "counts",
scale = "log") + ggtitle("Compensated")
plot1 + plot2
The following transformations are possible: arcsinh and biexp.
Transformation is performed on “counts” data of the DefaultAssay of the Seurat object. The DefaultAssay is either “fcs” or “comp” depending on whether compensation was performed or not.
The transform_data function will return a Seurat object with transformed data written into the “data” slot.
Both arcsinh and biexp transformation can and should be used with custom parameters. These parameters should be provided in the panel data.frame, which is stored in seu@misc upon Seurat object creation.
seu <- transform_data(seu, "arcsinh")
plot1 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "fcs", slot = "counts", scale = "log") +
ggtitle("Untransformed")
plot2 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", slot = "data") + ggtitle("Arcsinh transformation")
plot1 + plot2
Flow- and Mass cytometry data can contain millions of cells. Downsampling is optional but should be done at this step of the pipeline if desired.
Depending on computing power, it is advisable to downsample data. Furthermore, if number of cells differs largely between samples, it is advisable to downsample so overall differences are not driven by individual samples with high number of cells.
It is recommended to save the original non-downsampled object for potential applications later.
# set seed for downsampling for reproducibility
set.seed(42)
# downsample using Seurats subset function
seu <- subset(seu, downsample = 20000)
seu
## An object of class Seurat
## 89 features across 198659 samples within 3 assays
## Active assay: comp (39 features, 0 variable features)
## 2 other assays present: fcs, unused
# plot cellnumbers per sample
plot_cellnumber(seu)
Cytometry is an inherently noisy technology. There will be events close to the axis or artefacts due to antibody aggregation. Seumetry provides multiple quality control tools to achieve clean data for downstream analysis.
It can occur with cytometry, that some events have extremely negative or positive values. Here, we can remove these outlier events 1) by setting a manual threshold for each channel or 2) by using an automatic removal algorithm based on isolation forest.
Use a named vector to indicate threshold for each channel. The function can take positive or negative thresholds.
# check each channel to determine if they contain outlier events
plot_cyto(seu, "CD4", "CD3", style = "2d_density", assay = "comp", slot = "data")
# manually set thresholds based on cyto plots
thresholds_neg = c(CD3 = -4, CD4 = -2)
thresholds_pos = c(CD3 = 4, CD4 = 4.5)
# remove values above or below thresholds
seu_manual = remove_outliers_manual(seu, thresholds_neg)
seu_manual = remove_outliers_manual(seu_manual, thresholds_pos, negative = FALSE)
plot_cyto(seu_manual, "CD4", "CD3", style = "2d_density", assay = "comp", slot = "data")
Outliers are detected using an isolation forest. In ungated flow cytometry data, this algorithm will mainly remove axis-near events, dead cells, and doublets.
The threshold at which score an event is regarded an outlier can be supplied. The threshold is usually between 0.6 and 0.8. Default threshold: 0.7.
The function will save two features for each cell into meta.data of the Seurat Object: 1) the isolation forest score (seu\(outlier_score) and whether an event passed the threshold or not (seu\)outlier).
# identify outliers (default threshold: 0.7)
seu <- detect_outliers(seu)
plot1 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", slot = "data") + ggtitle("Before removal")
plot2 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", slot = "data", style = "point",
color = "outlier_score") + ggtitle("Outlier score")
plot3 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", slot = "data", style = "point",
color = "outlier") + ggtitle("Outlier")
plot4 <- plot_cyto(subset(seu, subset = outlier == FALSE), x = "CD4", y = "CD3",
assay = "comp", slot = "data") + ggtitle("After removal")
plot1 + plot2 + plot3 + plot4
Remove outliers from Seurat object.
# remove outliers from Seurat object
seu <- subset(seu, subset = outlier == FALSE)
With high-dimensional Flow cytometry, antibody aggregates can occur. These are usually characterized by highly co-linear events that form diagonal structures in XY plots. This algorithm can identify potentially problematic channel combinations and identify aggregates in these channels.
Here is an example of such aggregates:
plot1 <- plot_cyto(seu, x = "CD28", y = "CXCR5")
plot2 <- plot_cyto(seu, x = "IL15-2RB", y = "CD123")
plot1 + plot2
The first step is to identify channel combinations that potentially contain aggregates. This is done by using only double positive events (default: events > 1) and running a pearson correlation.
The function “detect_aggregate_channels” returns a list of 2 matrices and 1 data.frame.
Matrix 1: pearson correlation matrix.
Matrix 2: binary correlation matrix based on threshold of pearson R (default threshold = 0.7).
Data.frame contains the channel combinations that passed the pearson R threshold.
problem_channels <- detect_aggregate_channels(seu)
These correlation matrices can be plotted, for example, using a heatmap.
pheatmap::pheatmap(problem_channels[[1]], main = "Pearson R")
pheatmap::pheatmap(problem_channels[[2]], main = "Pearson R > 0.7")
Next, we can plot the channel combinations that potentially contain aggregates and manually assess if these channels are really problematic.
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
"Channel_1"], y = problem_channels[[3]][i, "Channel_2"])
do.call(gridExtra::grid.arrange, c(plots, ncol = 3))
If a channel combination does not look like it contains aggregates, remove that row from the data.frame (problem_channels[[3]]).
Next, we can detect aggregates in the problematic channels using a modified RANSAC algorithm. Each cell will get an aggregate_score, which is the number of channel combinations in which it was labelled as an aggregate. By default, potential aggregates are only labelled real aggregates if they occur in at least 2 channel combinations. The aggregate_score and whether a cell has passed the aggregate_score threshold (default: 2) is stored in meta.data of the Seurat object.
seu <- detect_aggregates(seu, problem_channels[[3]])
# Check the aggregate removal performance
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
"Channel_1"], y = problem_channels[[3]][i, "Channel_2"], style = "point", color = "aggregate_score")
do.call(gridExtra::grid.arrange, c(plots, ncol = 3))
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
"Channel_1"], y = problem_channels[[3]][i, "Channel_2"], style = "point", color = "aggregate")
do.call(gridExtra::grid.arrange, c(plots, ncol = 3))
If the result is satisfactory, the aggregates can be removed.
seu <- subset(seu, subset = aggregate == FALSE)
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
"Channel_1"], y = problem_channels[[3]][i, "Channel_2"])
do.call(gridExtra::grid.arrange, c(plots, ncol = 3))
It is possible that some aggregates still remain. In that case it is advisable to fine tune the parameters of the aggregate_channels function.
Alternatively, remaining artefacts may be removed by detecting outliers again using the isolation forest.
seu <- detect_outliers(seu, score_threshold = 0.6)
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
"Channel_1"], y = problem_channels[[3]][i, "Channel_2"], style = "point", color = "outlier")
do.call(gridExtra::grid.arrange, c(plots, ncol = 3))
# remove outliers from Seurat object
seu <- subset(seu, subset = outlier == FALSE)
UMAP, clustering, plotting etc is completely compatible with Seurat workflow. If UMAPs and clustering do not separate celltypes well, it can help to only select features for PCA that distinguish expected cell subsets, e.g. CD3, CD4, and CD8 if working with T cells.
Here, all features are used for scaling, PCA. For UMAP and clustering, all PCs are used.
# scale data
seu <- ScaleData(seu, features = row.names(seu))
# run PCA for all features
seu <- RunPCA(seu, features = row.names(seu), approx = FALSE)
# run UMAP using all PCs
seu <- RunUMAP(seu, dims = 1:length(seu@reductions$pca))
# run clustering using all PCs
seu <- FindNeighbors(seu, dims = 1:length(seu@reductions$pca))
seu <- FindClusters(seu, resolution = 0.1)
Now, all Seurat visualisations and downstream processing tools are available. For more detail, see https://satijalab.org/seurat/
plot1 <- DimPlot(seu, label = TRUE) + set_theme
plot2 <- DimPlot(seu, group.by = "sample_id") + set_theme
plot3 <- DimPlot(seu, group.by = "layer") + set_theme
plot1 + plot2 + plot3
VlnPlot(seu, features = c("CD8", "CD4"), pt.size = 0) & set_theme
It is recommended to use a custom colorpalette when using UMAP to plot expression of markers.
FeaturePlot(seu, features = c("CD8", "CD4")) & scale_color_gradientn(colors = gradient_colors) &
set_theme
Not only visualisations, but also other Seurat features are fully functional. For example, FindAllMarkers can be used to help annotate clusters.
markers <- FindAllMarkers(seu, only.pos = TRUE)
top5 = markers %>%
group_by(cluster) %>%
top_n(n = 5, wt = avg_log2FC)
DotPlot(seu, unique(top5$gene), assay = "comp") + scale_color_gradientn(colors = gradient_colors) +
theme_linedraw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
RotatedAxis()
Differential abundance (DA) analysis is implemented in the Seumetry workflow based on the edgeR package. DA is done using a generalized linear model (GLM) to model number of cells per group given attribute/feature x condition (e.g. seurat_clusters x treatment) and a likelihood ratio rest (LRT) to make a contrast between all conditions.
Using the function with “check_coef = TRUE” will return all coefficients to design the proper contrast for the DA analysis.
# first check the contrasts
differential_abundance(seu, attribute = "seurat_clusters", group_by = "sample_id",
formula = as.formula("~0+layer"), check_coef = TRUE)
## [1] "layerIEL" "layerLPL"
Here, we compare IEL vs LPL using a contrast of c(1, -1).
# calculate differential abundance: UC vs healthy
da_res <- differential_abundance(seu, attribute = "seurat_clusters", group_by = "sample_id",
formula = as.formula("~0+layer"), contrast = c(1, -1))
head(da_res)
## logFC logCPM LR PValue FDR
## 10 -2.6080289 11.53333 9.1669816 0.002464207 0.02609839
## 4 -2.7789907 14.69200 7.9740892 0.004745163 0.02609839
## 3 -3.1508768 15.63377 3.8023178 0.051181688 0.18766619
## 8 -1.5108192 13.63196 2.9062036 0.088239359 0.24265824
## 6 -2.3122751 13.91885 1.3506560 0.245163459 0.53935961
## 0 0.4649679 19.05763 0.5806879 0.446042812 0.74285370
The results can be exported as a table and/or plotted, for example, using EnhancedVolcano.
EnhancedVolcano::EnhancedVolcano(da_res, lab = rownames(da_res), drawConnectors = TRUE,
x = "logFC", y = "FDR", title = "IEL vs LPL", pCutoff = 0.05, FCcutoff = 0.2,
ylim = c(0, 2.5), xlim = c(-4, 4)) + set_theme
In addition to Seurat “FindMarkers” function, Seumetry also includes a feature to find differential expression (DE) of markers based on a pseudobulk analysis.
To this end, the median fluorescent intensity is calculated using the Seumetry function “median_expression”.
The DE analysis is based on the limma package and a linear model is fit followed by empirical Bayes statistics for differential expression (eBayes).
de_res <- DE_pseudobulk(seu, fixed_vars = "layer", contrast = "layerIEL-layerLPL")
head(de_res)
## logFC AveExpr t P.Value adj.P.Val B
## CCR6 0.12147671 0.0534059 2.821960 0.01915969 0.7472279 -3.160006
## HLA-DR -0.19894703 0.2026373 -2.035465 0.07094009 0.9593021 -3.939721
## CXCR3 0.35726377 0.6875903 1.934553 0.08364929 0.9593021 -4.039075
## CD8 0.82202071 1.9083093 1.568040 0.14985778 0.9593021 -4.387271
## CD103 0.42022311 1.9011594 1.309221 0.22153364 0.9593021 -4.613302
## PD1 -0.05660534 0.1746178 -1.299519 0.22470753 0.9593021 -4.621335
The results can be exported as a table and/or plotted, for example, using EnhancedVolcano.
EnhancedVolcano::EnhancedVolcano(de_res, lab = rownames(de_res), drawConnectors = TRUE,
x = "logFC", y = "P.Value", title = "IEL vs LPL", pCutoff = 0.1, FCcutoff = 0.2,
xlim = c(-1, 1), ylim = c(0, 2)) + set_theme
Sometimes it is useful to create cytometry-like plots such as density plots. Seumetry includes various plotting functionalities using “plot_cyto”.
plot_cyto(seu, x = "CD3", y = "CD4")
plot_cyto(seu, x = "CD3", y = "CD4", style = "point", color = "sample_id")
plot_cyto(seu, x = "CD3", style = "density", color = "sample_id")
A principal component analysis (PCA) based on median marker expression per sample can give a good impression of differences between conditions or replicates. Seumetry includes the “plot_pca” function to compute a sample-level PCA. The PCA can also be grouped by other factors than sample (based on seu@meta.data).
plot_pca(seu)
plot_pca(seu, group_by = "sample_id", color = "layer")
Seumetry includes a function to visualise the frequency or proportion of cells based on metadata columns. For example, it can be useful to assess whether the frequency of clusters is stable across different samples or differents between conditions.
plot_frequency(seu, "seurat_clusters", "sample_id")
plot_frequency(seu, "seurat_clusters", "layer")
plot_cellnumber(seu, "sample_id")
The Seumetry function “median_expression” can be used similar to Seurat “average_expression” function, but uses the median instead of the average. This function can be useful, for example, to plot median expression per cluster or condition.
sample_mfi <- median_expression(seu, group_by = "sample_id")
head(sample_mfi)
## 303_Ileum_IEL 303_Ileum_LPL 304_Ileum_IEL 304_Ileum_LPL 306_Ileum_IEL
## IL15RA -0.17426727 -0.20895809 0.08866171 0.11169343 -0.067285577
## CD27 -0.02206521 -0.05455738 0.01819752 -0.03724345 0.005262813
## CCR7 -0.06661292 -0.12887841 -0.03415701 -0.04021653 -0.041587133
## CD127 0.88077292 1.30208107 0.61359709 0.87511319 0.626864898
## CD1C 0.09538533 0.12469144 -0.18241319 -0.18354583 -0.046564421
## CD141 -0.28070071 -0.68942647 -0.05885200 -0.13515761 0.117691353
## 306_Ileum_LPL 309_Ileum_IEL 309_Ileum_LPL 310_Ileum_IEL 310_Ileum_LPL
## IL15RA -0.09089852 -0.10834529 -0.171781019 -0.10533218 -0.15235189
## CD27 -0.01122545 0.02059490 0.002202564 0.01028800 0.10182264
## CCR7 -0.03377532 -0.04077703 -0.049045663 -0.08163594 -0.02346143
## CD127 0.73624058 0.63497176 0.674176444 0.45260544 0.17213334
## CD1C -0.03263402 0.01600431 0.095684599 0.02203863 0.26617907
## CD141 0.15456455 -0.23181249 -0.248818314 -0.64322970 -0.52487300
Sometimes it can be useful to view specific clusters in external cytometry Software. For this purpose, Seumetry includes a function to export FCS files. This function is based on the initial “fcs_fs” object created for loading FCS files.
# subset Seurat object to specific cluster
seu_sub <- subset(seu, subset = seurat_clusters == 1)
# export FCS file containing cells from this cluster
export_fcs(fcs_fs, cell_ids = colnames(seu_sub), filename = "cluster_1.fcs")